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^ I Abstract 

J> I The spin-dependent corrections to the static inter-quark potential are phenomeno- 

■ logically relevant to describing the fine and hyperfine spin splitting of the heavy 
. quarkonium spectra. We investigate these corrections, which are represented as the 

0\ ] field strength correlators on the quark-antiquark source, in SU(3) lattice gauge the- 

' ory. We use the Polyakov loop correlation function as the quark-antiquark source, 

. and by employing the multi-level algorithm, we obtain remarkably clean signals for 

these corrections up to intermediate distances of around 0.6 fm. Our observation 
suggests several new features of the corrections. 



1 Introduction 



X ■ The spin-dependent potentials are parts of relativistic corrections to the 

■ static quark-antiquark potential, which depend on quark spin, and are phe- 

nomenologically relevant to describing the fine and hyperfine splitting of heavy 
quarkonium spectra [T1I2I31H] . Thus it is interesting to address these corrections 
from QCD and to compare with the observed spectra. 

The relativistic corrections are usually classified in powers of the inverse of 
heavy quark mass m (or quark velocity v) and it is well-known that in QCD the 
leading spin-dependent corrections show up at 0(l/m^) [51I61I71I81I91IT0] . These 
spin-dependent corrections were also derived systematically within an effective 
field theory framework called potential nonrelativistic QCD (pNRQCD) [llL 
pNRQCD is obtained by integrating out the scales above m S> Aqcd in QCD[^ 
first, which leads to NRQCD and then mv, leaving a typical scale of 

the binding energy of heavy quarkonium mv"^ [T¥lll5fl6] . 

The spin-dependent potential is summarized in the form 



^QCD is assumed to be a few hundred of MeV 
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+ (c?cgV4(r) - A8nCpasdJ^'\r)) , (1.1) 

where fi and r2 = l^i ~ "'^2!) denote the positions of quark and anti- 
quark, uii and 7712 the masses, si and S2 the spins (s = a/2 with a being 
the Pauh matrices), and /i = —I2 = I the orbital angular momenta. Vo(r) 
is the spin-independent static potential at 0{m^) and the prime denotes the 
derivative with respect to r. V{{r), V2{r), Vsi^r) and V4(r) are the functions 
which depend only on r. In what follows we call these functions loosely the 
spin-dependent potentials. Cp{iJ,,mi) {i = 1,2) is the matching coefficient in 
the (p)NRQCD Lagrangian which multiplies the term a ■ B/{2mi) and this 
coefficient plays an important role when connecting QCD at a scale /i with 
(p)NRQCD at scales rrii. We have defined {c^f'^ ±cP)/2. For equal 

quark and antiquark masses {mi = 7772), cf vanishes ^J. When 

the matching is performed at tree-level of perturbation theory, the coeffi- 
cient is c'p = 1 [17] and then Eq. (11. ip is reduced to the expression given 
in Refs. [5]|6f8] . as = j (47r) is the strong coupling and Cp = 4/3 the Casimir 
charge of the fundamental representation, and the mixing coefficient of the 
four-quark operator in the (p)NRQCD Lagrangian (see e.g. the Appendix E 
of Ref. 0). 

Given the field strength Fi,^, where the color-electric and the color- 
magnetic fields are defined by Ei = F^i and Bi = eijkFjk/2, respectively[fj 
the spin-dependent potentials in Eq. (11.11) are expressed as 

'^V;ir)=e,,k\m^J^ dtt{{g^B,{n,ti)E^{n,t2))) , (1.2) 
7^2(0 = e^jfc Jim ^"dt t{{g^B,in,ti)Ej{f2,t2))) , (1.3) 



^ Throughout this paper we work in Euclidean space and assume that the repeated 
spinor (Latin) and color (Greek) indices are summed over from 1 to 4 and from 1 
to 3, respectively, unless explicitly stated. 
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- f )^3(r) + fV,{r) = 2 Jim J^\t {{g' B,{nM)B,{r,M))) {I A) 

Here t = t2 — ti denotes the relative temporal distance between two field 
strength operators. The double bracket ((■■■)) represents the expectation value 
of the field strength correlator, where the field strength operators are attached 
to the quark- ant iquark source in a gauge invariant way. In Refs. |5|l6|l8] . these 
expressions were given in the double-integral form with the Wilson loop, which 
can be reduced to the single-integral form through the spectral representation 
of the field strength correlators derived from the transfer matrix theory. How- 
ever, it should be noted that the authors of Ref. [11] pointed out that one of 
the spin-orbit potentials ^'(r) in Refs. [5|6|8] was underestimated by a factor 
two. The expressions in Eqs. fll.2l) - fll.4p are consistent with Ref. [H] apart 
from the space-time metric; here we employ the Euclidean metric, while the 
Minkowski metric is used in Ref. [11] . Ll| 

As the expressions of the spin-dependent potentials in Eqs. f ll.2p - fll.4l) 
are essentially nonperturbative, these potentials can be studied in a frame- 
work beyond perturbation theory, for instance, by using lattice QCD Monte 
Carlo simulations. Rather, as the typical scale of the momentum mv can be of 
the order as Aqcd due to nonrelativistic nature of the system t> ^ 1, it is not 
clear a priori that the perturbative determination of the potential is justified, 
and indeed, nonperturbative contributions are expected in the spin-orbit po- 
tentials Vl{r) and V2(r); they are related to the static potential through the 
Gromes relation [JHUH], i-e. V^ij-) = V"2(r) — V({r), where Vo(r) is known to 
contain a nonperturbative long-ranged component characterized by the string 
tension. This relation was derived by exploiting the Lorentz invariance of the 
field strength correlators, which does not depend on the order of perturbation 
theory. 

The determination of the spin-dependent potentials using lattice QCD 
simulations goes back to the 1980s [2nll2T]l22]l23ll2if25]^^ and to the 
1990s [28]l29f30j . The qualitative findings (quantitative to some extent) in 
these earlier studies indicated that the spin-orbit potential V({r) contains the 
long-ranged nonperturbative component, while all other potentials are relevant 
only to short-range physics as expected from the one-gluon exchange interac- 
tion. However, one observes that the spin-dependent potentials (in particular, 
the spin-orbit potential V{) of even the latest studies [29|5D] suffer from large 
numerical errors, which can obscure their behavior already at intermediate 
distances. For phenomenological applications of these potentials, it is clearly 
important to determine their functional form as accurately as possible. 

In the present paper we thus revisit the determination of the spin- 
dependent potentials with lattice QCD within the quenched approximation, 
aiming at reducing the numerical errors with a new approach. There are mainly 



^ The change of metric from Minkowski to Euclidean space-time is achieved by 
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two possible sources of numerical errors apart from the systematic error due 
to discretization of space-time. One is the statistical error for the expectation 
value of the field strength correlator, and the other is the systematic error 
associated with the integration over r and the extrapolation of r ^ oo in 
Eqs. f ll.2l) - fll.4p . In order to control the total error, first of all, one needs to 
evaluate the field strength correlator precisely, as otherwise its uncertainty is 
enhanced in the following procedures. 

Our idea is then to employ the multi-level algorithm [3T|32] for measuring 
the field strength correlators |33fM] , as we expect clean signals even at larger r 
and t. This algorithm also allows us to use the Polyakov loop correlation 
function (PLCF), a pair of Polyakov loops P separated by a distance r, as the 
quark-antiquark source instead of the Wilson loop. In fact, if one relies on the 
commonly employed simulation algorithms, it is almost impossible to evaluate 
the field strength correlators on the PLCF, or the PLCF itself, at intermediate 
distances with reasonable computational effort, since the expectation value of 
the PLCF at zero temperature is smaller by several orders of magnitude than 
that of the Wilson loops and is easily hidden in the statistical noise. However, 
as we will show in the next section, if one manages to obtain accurate data 
for the field strength correlators on the PLCF, the systematic errors from the 
integration and the extrapolation can be avoided. The key idea is to employ the 
spectral representation of the field strength correlators on the PLCF, which 
is plugged into Eqs. (fL2!) - f|Li|) . 

This paper is organized as follows. In sect. [2], we describe our procedures 
to obtain the spin-dependent potentials, which contain the derivation of the 
spectral representation of the field strength correlators on the PLCF and of 
the spin-dependent potentials, the definition of the field strength operators on 
the lattice, and the implementation of the multi- level algorithm. In sect. [3l we 
show numerical results, followed by analyses and discussions. The summary 
is given in sect. HI In this paper, we will not discuss the matching coefficient 
but the interested reader can refer to the discussion in Ref. P]. We also plan 
to revisit this issue in our future studies. 



2 Numerical procedures 

In this section, we describe the spectral representation of the field strength 
correlators on the PLCF and of the spin-dependent potentials. We provide 
the definition of the field strength operators on the lattice, and explain the 
implementation of the multi-level algorithm. The standard Wilson action is 
most preferable for this algorithm because its action density is locally defined 
by plaquette and thus we shall use this action in our present simulation. The 
lattice volume is L^T and periodic boundary conditions are imposed in all 
directions. 
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2.1 Spectral representation of the field strength correlators and of the spin- 
dependent potentials 



Let us derive the spectral representation of the field strength correlators 
on the PLCF using the transfer matrix formalism. We follow the notation used 
in Ref . [32] , in which the spectral representation of the PLCF is discussed. We 



consider the transfer matrix in the temporal gauge T = g-™ ^j-^j^j-^ a^^g on 
the states on the space of all spatial lattice gauge fields at a given time, 
where a denotes the lattice spacing. We also introduce the projectors P onto 
the subspace of gauge- invariant states and P3(g)3*(ri, to the subspace of the 
states in the 3 ® 3* representation of SU(3). Then the partition functions in 
the sector corresponding to P and P3(g)3*(ri, r2) are given hy Z = TrlPe^"^"^} 
and 2^3®3*(n,^2) = |Tr{P3^3.(fi,f2)e"'^^}, respectively. 

Firstly, we consider the spectral representation of a double-bracket corre- 
lator for operators Oi(ti) and 02(^2), which are attached to either side of the 
PLCF (the same side or the opposite side). 



((Oi(ti)02(t2))) 



(Oi(ti)02(t2))p(ri)P*(r2) 



(P(fi)P*(f2)) 
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p. 



Z 

303*(n,r2)e" 



-23^3* (n, ^2) 



(2.1; 



where we have used the identity {P{fi)P*{r2)) = Z303. (ri, r2)/Z. Inserting 
the complete set of eigenstates in the 3 ® 3* representation \n) = \n;ri,r2), 
which satisfy T\n) = e~^"^^^"'\n) with energies En{r) > 0, we obtain 



((Oi(ti)02(t2))) 



E^o,m=o {n\0i\m){m\02\n) e-^^'e--^"^^-*) 



^n=0 ^ 



-EnT 



(2.2) 



We denote the energy gap between two eigenstates as AEmn{r) = Emir) — 
En{r). Then, up to terms involving exponential factors equal to or smaller 
than e-(^-^io)'^, Eq. (Q is reduced to 

((Oi(ti)O2(t2))) = (0|Oi|0)(0|O2|0) 

+ E ( (0|Oi|m)(m|O2|0)e-('^^'"«)* + (m|Oi|0)(0|02|m)e-(^^-o)(^-*M 



m=l 



+0(e 



(2.3) 



In the case of the field strength correlators, we can further simplify 
Eq. (12.31) by using the properties of the color-magnetic and color-electric field 
operators under the time reversal; we have relations 
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(ig^B,{t,)E,{t,))) = -{{g'B,{h)E,{h))) , (2.4) 
{{9'B,it,)B,it,))) = {{g'B,{t,)B,{t,))) , (2.5) 

which, for the matrix elements, read 

{m\gB,\Q){Q\gE^\m) = -{Q\gB,\m){m\gE^\Q) , (2.6) 
{m\gB,\Q){Q\gB^\m) = {Q\gB,\m){m\gB,\Q) , (2.7) 

for m > 1. Moreover, {Q\gBi\Q) = since Bi is odd under CP transformations. 
The field strength correlators in Eqs. fll.2l) - fll.4p are thus expressed as 

oo 

{{g^B,{futi)E,in,t2)))=2Y.{0\gBiin)\m){m\gE,in)\0) 



m=l 

^^-iAE^o)T/2 sinh((AE„o)(r/2 - 1)) 

^ -(A£;lo)T^| 

oo 



^Q(g-(AMo)i) ^ (2. 

oo 

{{g^B,{n,t^)Ejif2,t2)))=2Y,{0\gBiin)\m){m\gE,ir2m 



m=l 

xe-^^^"'°^^/^smh{{AE^o){T/2 - 1)) 

^ -(Ai^lo)T^| 

oo 



+0(e-^^^^«^^ ) , (2.9) 

oo 

((/5,(fi,ti)5,(f2,t2)))=2 5:(0|(75,(fi)|m)(m|(75,(f2)|0) 

m=l 

^g-(A£;„o)T/2 cosh((A£;„o)(T/2 - 1)) 
+0(e-(^^^«)^) . (2.10) 

After inserting these expressions into Eqs. fll.2l) - fll.4p . we can carry out the 
integration and extrapolation, which imply that 

Mm dt---= Mm dt--- (2.11) 

with an arbitrary rj G (0, 1/2]. Thereby we obtain the spectral representation 
of the spin-dependent potentials, which consists of the matrix elements and 
the energy gaps. 

For the simplest parametrization fi = = (0, 0, 0) with ti = and 
r2 = r = (r, 0, 0) with t2 = t, which is the actual setting of our simulation, we 
have explicitly 

y,^^^ ^ 2 £ i^\9By{0)\m){m\gEM\0) ^ ^2.12) 
^ (0|^i?,(0)|m)(m|^E,(r)|0) 

m=l l^-C/mOj 
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{Q\gB,{Q)\m){m\gB,{f)\{)) {Q\gBy{Q)\m){m\gBy{r)\Q) ' 

(2.14) 

(0|^i?.(0)|m)(m|^i?.(r)|0) ^ ^ (0|^j?,(0)|m)(m|^j?,(f)|0) l 

(2.15) 

where we have used the relations 

((^72fi^(0,0)E,(0,t))) = -((/5,(0,0)E,(0,t))) , (2.16) 
((/5,(0,0)E,(f,t))) = -((/i?,(0,0)E,(f,t))) , (2.17) 
((/i?,(0,0)S,(f,t))) = ((/i?,(0,0)i?,(f,t))) . (2.18) 

We note that the error term in the field strength correlator of 0(e~^^'^^°''^) in 
Eqs. ( I2.8p - fl2.10l) is assumed to be negligible, which is the case for large T. 

Now our procedure to compute the spin-dependent potentials is as follows; 
we evaluate the field strength correlators for various r and t, fit them to the 
spectral representation in Eqs. fl2.8p - fl2.10p . thereby determining the matrix 
elements and the energy gaps, and insert them into Eqs. f l2.12p -f l2.15p . Here, 
the hyperbolic sine or cosine function in Eqs. fl2.8p -f l2TT0l) is typical for the 
PLCF, which automatically takes into account the effect of the finite temporal 
lattice size in the fit. 

Note that if one uses the Wilson loop at this point, the spectral repre- 
sentation is just a multi-exponential function and the leading error term is of 
QI^^-{AEw){^t)-^^ where At is the relative temporal distance between the spatial 
part of the Wilson loop and the field strength operator [29.30J. Denoting the 
temporal extent of the Wilson loop by T^, one can fit the data in the range 
t G [Q,Tw — 2 At], where is at most T/2 because of the periodicity of the 
lattice volume. Clearly the available fit range is more restricted than in the 
PLCF case, even if At/a is chosen as small as possible, say one or two. It may 
be possible to suppress the error term by applying smearing techniques to the 
spatial links. However, it is not immediately clear if this procedure really cures 
the error term. At least, one needs fine tuning of the smearing parameters and 
further systematic checks. 

2.2 Field strength operator on the lattice 

On the lattice, we use the field strength operator defined by ga^F^y{s) = 
(U^uis) - Ul,{s))/{2i), where U^,{s) = U,{s)U,{s + /i)[/t(s + z>)f/t(s) is the 
plaquette variable at a site s = {sg, St) with a spatial site Sg and a temporal 
site St. We also define U_fj_{s) = Uj^{s — fi). Practically, we construct the 
color-electric and color-magnetic field operators by averaging the field strength 
operator as 



V,{r) = 2 E 

m=l 



V,{r) = 2 Y: 

m=l 
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ga^Eiis) = ^ga' {F^^{s)+F_,^{s)) , (2.19) 
ga^B^is) = \gahi,k {Fjk{s) + Fk-j{s) + + , (2.20) 

o 

where we assume that Ei{s) is defined on (s^, St + 1/2), and Bi{s) on (ss, Sj), 
respectively. 

Now, as seen from Eqs. fl2.12p - fl2.15p . the spin-dependent potentials con- 
sist not only of the energy gap but also of the matrix element of the field 
strength operator, and thus one needs to take into account the renormaliza- 
tion of the latter. This is due to the fact that the field strength operators 
depend explicitly on the lattice cutoff a. In the absence of a viable nonper- 
turbative renormalization prescription for the field strength operators in the 
presence of the quark-antiquark source, we follow here the Huntley-Michael 
(HM) procedure [23], which was also used in Refs. [29|30] . This procedure is 
inspired by the weak coupling expansion of the Wilson loop and is aimed at 
removing the self-energy contribution, at least, at 0{g'^). We define Ei and Bi 
from F^u{s) = {U^y{s) + f/^^(s))/2, and, by taking the average according to 
Eqs. f l^rg]) and ([22DD, we compute 

ZEXr) = l/{m. ZsXr) = l/m), (2.21) 

where E or B are attached to either side of the PLCF. These factors are then 
multiplied to the field strength operators in Eqs. fl2.19p and fl2.20p accordingly. 
Note that Ze, and Zb, determined in this way can depend on r and also on 
the relative orientation of the field strength operator to the quark-antiquark 
axis. 

One may find that the HM procedure is quite similar to tadpole improve- 
ment [33], where the corresponding renormalization factor is defined by the 
inverse of the expectation value of the plaquette variables, Ztad = '^/{Un), 
where 

which was used e.g. in Refs. [201128]. Indeed, if the factorization of the correlator 
{Ff_iu)pp* = {Ff_iu) {P P*) holdsj^ Zsi and Ze^ are reduced to ^tad- However, 
as was pointed out in [23], the tadpole factor does not remove the self energy 
even to 0{g^) if the correlator involves the electric field operator. 

2.3 Multi-level algorithm for the field strength correlator 

We now describe the multi-level algorithm [3T|32] for computing the field 
strength correlators, restricting the discussion to the lowest level. The essence 



Numerically, this factorization is approximately satisfied. 
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Table 1 

A minimal set of sublattice correlators for the static potential and the spin- 
dependent potentials. 



Potential 



Sublattice correlators 




V3,V^ 



of the multi-level algorithm is to construct the desired correlation function, 
which may have a very small expectation value, from the product of the rela- 
tively large "sublattice average" of its components. We will also refer to such 
a component as the sublattice correlator. 

The sublattice is defined by dividing the lattice volume into several layers 
along the time direction, and thus a sublattice consists of a certain number of 
time slices A^tsi (the number of sublattices is then A^sub = (T / a) / Ntsi, which 
is assumed to be an integer). The sublattice correlators are evaluated in each 
sublattice after updating the gauge field with a mixture of heatbath (HB) 
and over-relaxation (OR) steps, while the space-like links on the boundary 
between sublattices remain intact during the update. We refer to this pro- 
cedure as the "internal update". We repeat the internal update A^iupd times 
until we obtain a stable signal for the sublattice correlators. Next, we multi- 
ply these sublattice correlators in a suitable way to complete the correlation 
function, as described below. Thereby the correlation function is obtained for 
one configuration. We then update the whole set of links without specifying 
any layer to obtain another independent gauge configuration and repeat the 
above sublattice averaging. The computational cost for one configuration is 
rather large, but one can observe a signal already from a few configurations 
once iVtsi and iViupj are appropriately chosen. 

In the current simulation, the building blocks of the field strength corre- 
lators are the sublattice correlators listed in Table [H T represents the complex 
9x9 matrices that act on color tensors in the 3 ® 3* representation of SU(3). 
The subscripts of T in Table [1] denote the type of the sublattice correlators. 
The way of completing a sublattice correlator is as follows. Denoting the tem- 
poral sites as St = (^tsh'^sub), where itsi ^ runs within the extent of 
one sublattice labeled by isuh ^ [Ij^sub], a component of the Polyakov loop 
(timelike Wilson line V), the complex 3x3 matrices, in each sublattice is 
expressed as 



where we explicitly write the color labels in Greek letters. The direct product of 



( 



,itsi=i 



n Ut{Ss, isnhjtsl) 




(2.23) 
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two timelike Wilson lines V separated by a distance r is the simplest sublattice 
correlator, i.e. 

Tpp(Ss, isnh] r)af3^S = V{Ss, isnh)al3 ® V*{Ss + rX, isuh)-yS , (2.24) 

which is relevant to both the PLCF and the field strength correlators. 

Other sublattice correlators are constructed by inserting one or two field 
strength operators into the timelike Wilson line V. For instance, in order to 
obtain TpBy, "^pe^ etc., we compute in each sublattice the timelike Wilson 
line with a single insertion of the field strength operator and then its direct 
product with V. The argument of this type of correlators is (s^, ^sub; ^tsi), 
where itsi labels the timeslice where the field strength operator is inserted, 
its\ e [l,iVtsi]. The quantities Tp(^ByE,), ^ByE,, Tb^b, and TsySy represent the 
double-field-strength-operator-inserted sublattice correlators. The argument 
of these correlators is {ss,isuh',r,ired), where i^cd ^ [— A'tsi + 1, A^tsi — 1] is the 
relative temporal distance between two field strength operators. For Tp(^ByE^), 
two field strength operators are inserted into one of two timelike Wilson lines, 
while for TsyE^, Ts^p^ and TpyBy, they are inserted into both ones. 

The multiplication law of Tpp is then defined by 

{Tpp(s^, isub; r)Tpp{ss, isub + 1; r)}al3-yS 

= Tpp{Ss, isuh] r)ap'ya^Pp{Ss, ^sub + 1; t) pf^aS (2.25) 

and this multiplication law of color components is common to all other sub- 
lattice correlators. 

After taking the sublattice averages, we compute the PLCF for various r 

by 

PiO)P*ir) = J2 Re{[Tpp(s„ 1; r)] [Tpp(.„ 2; r)] x ■ • • 

x[Tpp(s,, A^sub - Ur)][Tpp{ss,Nsuh;r)]}aa-n ' (2-26) 

and the field strength correlators for various r and t by combining other sub- 
lattice correlators, where the translationally equivalent setting for space and 
time directions are averaged accordingly. Figure [T] illustrates the computa- 
tion of the field strength correlator {{g^By{ri,ti)Ez{ri,t2))) for V({r) (see also 
ref. [36] for a similar application of the multi-level algorithm, in which the 
electric-fiux profile between static charges was measured with the PLCF). 

In order to benefit from the multi-level algorithm, we need to optimize Atsi 
and A'iupd- They depend on the coupling /3 and on the distances to be investi- 
gated, which can be determined by looking at the behavior of the correlation 
function as a function of Ajupd for several iVtsi. As an empirical observation we 
note that aA^tsi = 0.3 — 0.4 fm is optimal in order to suppress the statistical 
fiuctuation of the correlation functions. 

In principle, one can apply the above computation to any direction of 
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T 



[ I t ] 



t 




Ntsl 



ri n 



r 



Fig. 1. How to complete {{g'^ By{ri,ti)Ez{ri,t2))) on the PLCF for V({r). Arrows 
at fi and r2 represent the Polyakov hnes for the static quark and antiquark. [• • • ] 
imphes the sublattice average. Other field strength correlators are obtained in a 
similar way. 



the quark- antiquark axis, r = (r^,., r^, r^,) with r = \Jr1 + ry + r'l. How- 

ever, one needs to take into account the fact that even with the simplest 
parametrization, r = (r, 0, 0), this algorithm requires a lot of memory to keep 
all T listed in Table [1] during the internal update. For a fixed distance and a 
fixed quark-antiquark axis using single precision, Tpp requires memory space 
{L/aY X A^sub X 162 x 4 bytes (= 1 work space unit wsu). Therefore, to com- 
pute all spin-dependent potentials in the same setting, one needs additionally 
SA^tsi wsu for the single- and 4(2A^tsi — 1) wsu for the double-field-strength- 
inserted sublattice correlators. For instance, on the 20^40 lattice with A^tsi = 4 
(^sub = 8), about 49 wsu = 2 GB memory is needed as 1 wsu = 42 MB. 

The way of obtaining the HM factors Ze^, and Zb^ using the multi- 
level algorithm is the same as above; we evaluate the sublattice average of 
Tp^^, Tp§^ and Tp^^ and complete the correlators in Eq. fl2.2ip . The addi- 
tional memory requirement is SiVtsi wsu in the above setting. 



3 Numerical results 

In this section, we present our numerical results. Simulation parameters 
are summarized in Table [2l We investigated the bare gauge couplings (3 = 6.0 
(a ~ 0.093 fm) on several lattice volumes and jS = 6.3 (a ~ 0.059 fm) on 
a 24^ lattice. The physical lattice spacing a was fixed in terms of the Sommer 
scale by setting tq = 0.5 fm [37]. One Monte Carlo update consisted of a 
combination of 1 HB and 5 OR steps. The number of internal updates iViupd 
for each (3 value was optimized empirically to obtain signals up to intermediate 
distances. We note that N^si and A^iupd could further be tuned so that even 
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Table 2 

Simulation parameters used in this study. The fourth column denotes the available 
quark-antiquark distances for the static potential Vq (before tree-level improve- 
ment), while [• • • ] applies to the spin-dependent potentials V(, V2, V3 and V4. 



(3 = 6/52 


a [fm] 


{L/af{T/a) 




r/a 






^iupd 


-^conf 


6.0 


0.093 


16^ 


2 - 


7 [2- 


6] 


4 


10000 


90 


6.0 


0.093 


20^ 


2 - 


9 [2 - 


6] 


4 


7000 


82 


6.0 


0.093 


20^40 


2 - 


9 [2 - 


7] 


4 


7000 


33 


6.3 


0.059 


244 


2 - 


8 [2 - 


6] 


6 


6000 


39 



larger distances become accessible. The statistical errors were estimated by 
applying the single elimination jackknife analysis. The various fit parameters 
were determined by minimizing defined with the full covariance matrix, 
and their errors were estimated from the distribution of the jackknife samples. 
For a consistency check we also evaluated the errors of the fit parameters 
from the minimum value of the through Ax^j^ = 1. In general, the errors 
from the jackknife analysis were the same or slightly larger compared to those 
estimated from Ax^j^^ = 1. 

3. 1 Static potential and its derivatives 

We first present the basic quantities extracted from the PLCF, i.e. the 
static potential and its derivatives with respect to the distance, in Figs. [2llU 



which are defined by 

Voin) = ln(P(0)P*(f)) + 0(e-(^^-)^) , (3.1) 

V^if) = ^{Voir)-Voir-a)} , (3.2) 

Ir'V^ir) = ^r^l {Vo{r + a) + K,(r - a) - 2K,(r)} = -c{r) . (3.3) 



We have applied tree-level improvement to the quark-antiquark distances in 
order to avoid an enhancement of lattice discretization effects especially at 
small distances pSll37f32j . so that the distances, r/, f and r are defined through 



the relations 

rfi = 47rG(r) , (3.4) 

= — {G{r - a) - G{r)} , (3.5) 
a 

f-3 = ^ lG(r + a) + G(r -a)- 2G(r)} , (3.6) 



where G{r) = G{r, 0, 0) is the Green function of the lattice Laplacian in three 
dimensions. For convenience we summarize these distances in Table RT] in the 
Appendix. 
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Fig. 2. Static potential Vo{ri) at (3 = 6.0 on the 20'^40 lattice and at /3 = 6.3 on 
the 24^ lattice. The constant term is subtracted. The dotted line is the fit curve 
Eq. p.7p . applied to the data at (3 = 6.0. 
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Fig. 3. The force VQ{r). The dotted line is the fit curve Eq. (j3.8p . applied to the 
data at /? = 6.0. 

We then fit the potential and the force (first derivative of the potential) 
to the functions 



c 

ar \- fi , 

r 
c 



(3.7) 
(3.8) 



and estimate the string tension a, the Coulombic coefficient c, and a con- 
stant II. The fit results at /? = 6.0 are summarized in Table [31 We find that 
the string tensions determined from either the potential or the force are con- 
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Fig. 4. c(f) = -f^V^'{r)/2. 



Table 3 

Fit results of the static potential and the force at /? = 6.0 on the 20^40 lattice with 
the fit functions in Eqs. (j3.7p and (j3.8p . The corresponding fit curves are plotted in 
Figs. El and El 





Fit range 


(r/a) 




c 


fia 


Xmin/^df 


^o(r) 


1.855 - i 


^.971 


0.0466(2) 


0.281(5) 


0.7169(5) 


3.8 




3.312 - i 


^.438 


0.0468(2) 


0.297(1) 




5.6 



sistent with each other. There is a small difference in c. However, this may be 
acceptable, since we see that c(f), which is extracted from the second deriva- 
tive of the static potential, is not strictly constant as a function of r as shown 
in Fig. m Thus c can be affected by the additional fit terms. Also, given that 
the estimates for c have not stabilized in the considered range of distances, it 
is not too surprising that a fit ansatz based on Eqs. (13.7!) and (13. 8p produces 
a large value of x^- Nevertheless, it is interesting to find that the fit curves 
characterize the global feature of the potential and the force. In these fits, 
we found no strong dependence on the fit range. We also examined the force 
obtained by the central derivative V^'(rc) = {Vo(r + a) — Vo(r — a)}/(2a), where 
fc is defined via = 47r{G(r — a) — G{r + a)}/(2a), and found the same 
curve as in Fig. [31 Later, we shall compare the values of a and c with those 
extracted from the spin-dependent potentials. 

At large enough distances, one may expect the value c = 7r/12 ^ 0.262 
from the bosonic string theory in four dimensions [32110] ■ However, as is clear 
from the plot of c(f) in Fig. HI we find c(f) 0.3 at r > 0.3 fm, which differs 
from this expectation by about 13 %. To accommodate a value of 7r/12 for 
the coefficient of the 1/r term, one would need higher order corrections at 
these distances. Note that the results for c(f) obtained here are consistent 
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Fig. 5. The HM factors at /3 = 6.0 (left) and (3 = 6.3 (right) as a function of r. 
The dashed hnes correspond to the tadpole estimate from the inverse of the expec- 
tation value of the plaquette: Ztad = 1-684 (([/□) = 0.59373(4)) at /3 = 6.0 and 
Ztad = 1.607 {{Un) = 0.62241(2)) at /3 = 6.3. 

with Refs. |32]HT] . 



3.2 HM renormalization factors 



In Fig. m we show the HM renormalization factors defined in Eq. (]2.2ip . 
together with the tadpole renormalization factor (see Table IA.2I in the Ap- 
pendix for the numerical values). Z^. are almost constant as a function of r, 
while exhibit some dependences on r and on the relative orientation of 
the operator to the quark- ant iquark axis at smaller distances. The HM fac- 
tors are generally smaller than the tadpole factor by less than 1 % for Z^. 
and at most 6 % for Z^. . Since the statistical fluctuations of these factors are 
much smaller than that of the field strength correlators, we ignore their errors 
when we multiply them to the field strength correlators, but take into account 
their r-dependences. Note that the observed tendencies are in agreement with 
Ref. |30], where the Wilson loop was used. 



3.3 Field strength correlators 



In Figs. [6] and [3, we show the various field strength correlators as a func- 
tion of t at /5 = 6.0 on the 20^ and 20^40 lattices, respectively, where r/a = 5 
is selected as an example. At smaller distances, the quality of the data is even 
better. Owing to the multi-level algorithm, the statistical accuracy of the data 
is unprecedented, which allows us see the typical behavior of correlators. 

We also include the fit curves in these figures, which are supplied by 
the spectral representation of the field strength correlators on the PLCF in 
subsection 12.11 We find that Eqs. fl2.8p - fl2.10p provide an excellent description 
of the behavior of the lattice data. Our fit procedure was as follows. We first 
folded the data of C{t) = {{g^By{6,0)E^{6,t))) and {{g^By{6,0)E^{f,t))) as 
{C{t)-C{T-t)}/2 C{t) with t e [0.5a, (T-a)/2], and the data oiC{t) = 
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{{g^B,{0,0)B,{f,t))) and {{g''By{0,0)By{r,t))) as {C(t) + C(T-t)}/2 ^ C{t) 
with t G [0,T/2]. Then, we fitted all available data points for each correlator 
in order to take into account as many excited states as possible in the spectral 
representations, except for t/a = 0.5 in {{g^By{0,0)E;;{0,t))), which was to 
avoid unwanted lattice effects due to the sharing of the same link variable in the 
two field strength operators. Here, since it is impossible to fix the parameters, 
matrix elements and the energy gaps, for all excited states, m > 1, with a finite 
number of data points, we truncated the expansion at a certain m = mmax- 
The validity of this truncation was verified by monitoring the values of and 
the integration results. We generally chose mmax such that Xmin/^df was of 
order 1, where A^^f is the number of degrees of freedom. All fit details and the 
fit results (i.e. the values of the integral) are tabulated in Tables IXl3] and PV.41 
in the Appendix. 

An important observation is that the integrals obtained for T = 20a and 
T = 40a are the same within errors, despite the fact that the behavior of 
the correlators around t = T/2 in Figs. [n]and[7|is obviously different. While 
the correlator computed for T = 20a is clearly distorted due to periodicity, 
this is not the case for T = 40a (note that the vertical axis is the same in 
both figures). This indicates that the hyperbolic sine or cosine function in the 
spectral representation provides a good description of the finite-T effect. In 
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other words, it is possible to extract the asymptotic value of the integral from 
the lattice with a relatively small value of T within this approach. At the 
same time, this confirms that the error term of 0{e~^^^^°^^) is negligible in 
this setting. We also examined a smaller lattice volume, 16^, at /3 = 6.0 and 
obtained the same result at r/a < 613 In this sense the spatial finite volume 
effect at r/a < 6 on the 20^ and 20^40 lattices is also negligible. 

The fit result at j3 = 6.3 on the 24^ lattice is listed in Table IA.5I in the 
Appendix. The finite volume effect at this /3 value is expected to be small, 
though we did not investigate this explicitly, since the physical size of the 
lattice volume is almost the same as for 16^ at /3 = 6.0. 

Finally, we point out several caveats, i) Although the fit works beautifully, 
one may not be able to assign a quantitative meaning to the resulting matrix 
elements and the energy gaps. This is because we truncated the spectral rep- 
resentation when performing the fit, and as a result, the lattice data, which in 
principle contains the contribution from all modes, are forced to be described 
by only a few modes. In this case, it is reasonable to regard only the value for 
the integral to be of quantitative significance, ii) As the energy gaps of the 
various field strength correlators in Eqs. fl2.8l) - fl2.10l) are the same, one may 
expect that the fit of each correlator provides the same energy gap, or one may 



These data are not presented in this paper but available on request. 



17 



attempt a simultaneous fit of all correlators with such a constraint. However, 
as the effective truncation level is not always common to all correlators, even 
at a fixed distance, which is also related to the fact that the matrix elements 
are not positive definite, this was not always the case, iii) The sensitivity of the 
fit result, namely the integration value, is mostly governed by the lowest en- 
ergy gap selected by the fit, which gives the dominant contribution at r ^ cxd. 
This is why we examined two lattice volumes with the same spatial size but 
different temporal extent and confirmed that such a systematic effect is neg- 
ligible. This fact supports our claim that the spectral representation of the 
field strength correlator is useful even though a truncation must necessarily 
be performed. 

3.4 Spin- dependent potentials 

The spin-dependent potentials, V{{r), ^'(r), V3(r) and V4(r) at P = 6.0 on 
the 20^40 lattice and at P = 6.3 on the 24^ lattice are presented in Figs. [HI El OS] 
and [131 respectively, expressed in physical units. These are the main results 
of this paper. Though we expect a scaling behavior for VsB{r) in Eq. (11. ip . 
both data at (3 = 6.0 and (3 = 6.3 for each potential seem to fall into one 
curve, which in turn may indicate that the matching coefficients should depend 
weakly on /3. The qualitative behavior of these potentials is not obscured by 
numerical errors. However, there is still room for improvement for the data 
with r > 0.3 fm at /5 = 6.3. The raw data in the lattice unit are summarized 
in Tables [Al6llA.8l in the Appendix^ The rest of this subsection is devoted to 
the interpretation of our data. In particular, we shall discuss the functional 
form of the dependence of the potentials on the distance r. 

We start by briefiy summarizing the theoretical expectation for the 
spin-dependent potentials. As mentioned in the introduction, Gromes de- 
rived a relation between the static potential and the spin-orbit potentials, 
VqIt) = V^'(r) — V({r), using the Lorentz (Poincare) invariance of the field 
strength correlators [T8|T9] . He also derived several inequalities for the spin- 
spin potentials based on refiection positivity, such as V^lr) > V4(r) and 
2Vs{r) + Vi{r) > [42]. These relations are nonperturbative, and can thus 
be checked directly on the latticejZ] Moreover, these relation do not depend 
on the matching scale. 

Another source of information comes from the systematic non-relativistic 
reduction of the Bethe-Salpeter (BS) equations within the instantaneous ap- 
proximation pp. Starting from the interaction kernel, which is assumed to be 

^ In these data, the HM factors are already multiplied, but the bare lattice data 
can be extracted by dividing the corresponding factors in Table IA.2I Starting from 
the bare data one can also test other renormalization procedures. 
^ One may of course expect a certain deviation from this relation on the lattice 
with a finite lattice cutoff a, since the strict Lorentz invariance is restored only in 
the continuum limit, a ^ 0. 
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Table 4 

The relation between the Lorentz structure of the effective kernel in the Bethe- 
Salpeter equation and the spin-dependent potentials [1]. 5(r), V{r) and P(r) are 
some scalar functions. If the interaction kernel has several components, the expected 
forms of the potentials are given by the sum of the corresponding terms. 



Kernel 


Voir) 


Vi{r) 


V2ir) 


Vsir) 


y^ir) 


Scalar 


S{r) 


-Sir) 











Vector 


V{r) 





Vir) 


-V"ir) +V'ir)/r 


2Ay(r) 


Pseudo-scalar 











P"ir) - P'(r)/r 


AP(r) 



a function of the norm squared of the relative momentum between a quark 
and an antiquark with various Lorentz structures, one arrives at a Breit-Fermi 
type effective Hamiltonian up to 0(l/m^) [^3|1^ . By comparing this effective 
Hamiltonian with Eq. (II. ip (where C'^p = 1 is assumed), one obtains the rela- 
tion between the Lorentz structure of the kernel and the spin-dependent poten- 
tials as summarized in Table HI This indicates that the Lorentz structure of the 
confining static potential can also be inferred from the structure of the spin- 
dependent potentials. For the special case of the one-gluon-exchange interac- 
tion, the kernel only consists of the Lorentz vector, and the spin-dependent 
potentials are explicitly given by 

^i'(r) = 0, Viir) = ^^, W = ^, F4(r) = 87rc5(3)(r), (3.9) 
where c = Cpag- 

We shall now investigate the r-dependence of our lattice data in more de- 
tail. We emphasize that, apart from the Gromes relation, no exact predictions 
exist for the behavior of the potentials beyond the short-distance regime. We 
have thus investigated the r-dependence of the potential by fitting the data to 
a particular model function, mostly guided by the short-distance predictions 
of Eq. (13. 9p . In cases where the latter clearly failed to describe the data, we 
have also resorted to effective parametrizations. The quality of a particular fit 
ansatz was judged by monitoring the value of Xmm/^df as computed using the 
full covariance matrix. Clearly, the underlying mechanism responsible for the 
observed behavior cannot be established rigorously in this manner. However, 
the main aim of this analysis is to provide nonperturbative input and guidance 
for future conceptual studies in this area. 

In the following we concentrate mostly on the dataset a,i (3 = 6.0, since 
it extends to larger distances compared with the data collected at /5 = 6.3. 
On the other hand, at smaller distances the results may be affected more 
by lattice artefacts. Indeed, around r 0.2 fm we occasionally observe small 
discrepancies for some potentials. Therefore, in the fits to the (3 = 6.0 dataset 
described below, we have mostly omitted the data point corresponding to the 
smallest separation r/a = 2. 
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Fig. 8. Spin-orbit potential V({r) at /? = 6.0 and (3 = 6.3. The dotted line is the fit 
curve Eq. (j3.10p . applied to the data of /? = 6.0. 
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Fig. 9. Spin-orbit potential V2(r) at /? = 6.0 and /? = 6.3. The dotted line is the fit 
curve Eq. (|3.8p . applied to the data of /3 = 6.0. 



We start our discussion with the spin-orbit potentials V({r) and V2(r). 
For V{{r), we find that the potential is negative and almost constant at r > 
0.25 fm (see Fig. [8]). This behavior clearly contradicts Eq. fl3.9p and suggests 
the existence of the Lorentz-scalar piece in the interaction kernel in terms 
of the BS equation. Our data at (3 = 6.3 suggest that one cannot exclude a 
deviation from a constant at small distances, an observation which was also 
made by Bali et al. [29!'30]. For V"2(r), we see a decreasing behavior with r (see 
Fig. [9]), which is not restricted to the short range, but rather seems to have a 
finite tail up to intermediate distances. 

Before discussing the functional form of V{{r) and V2{r) quantitatively, we 
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Fig. 10. Comparison between the force VJ)'(r) and the difference of the spin-orbit 
potentials V^{r) - y/(r) at /3 = 6.0 on the 20^40 lattice and (3 = 6.3 on the 24*^ 
lattice with the physical scale. The Gromes relation is V^'(r) = ^^'(r) — V({r). The 
dotted line is the fit curve for the force Vq, while the dashed line is for ~ ^i) 
where Eq. (j3.8p is used in both cases. 
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Fig. 11. The relative deviation from the Gromes relation, 1 — (V'2 — ^lOfit/^'-gt) 
(3 = 6.0. The fit functions are the same as in Fig. [10] and the dotted lines denote 
the 1(7 error band. 



may ask if these spin-orbit potentials satisfy the Gromes relation, since other- 
wise it is unclear whether any quantitative arguments make sense. In Fig. [TUl 
we compare the static force, VJ)'(r), with the difference of the spin-orbit po- 
tentials, V2(r) — V({r). We find quite a good agreement, indicating that the 
Gromes relation is apparently satisfied. We can examine this relation in more 
detail by fitting the difference V"2'(r) — V-[{r) at /3 = 6.0 to the same functional 
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form as the force in Eq. (13. 8p . Thereby we obtain a^2io^ = 0.0426(9) and 
c = 0.293(9) with XminZ-^df = 0.26. The string tension extracted in this way 
is about 8 % smaller than that in VqIt), while the Coulombic coefficients are 
in agreement within errors. We also plot the quantity 1 — (V2 — ^i)fit/K)fit 
Fig. [m From this we conclude that the Gromes relation is satisfied within 
(8 ± 1) % accuracy at r ^ 0.5 fm. Note that without the renormalization 
factor for the field strength operator, one would observe a strong deviation 
from the Gromes relation by a factor ^ 2.7 at f3 = 6.0. In this sense, the 
renormalization of the operator is crucial for satisfying the Gromes relation 
within a few percent level, especially when the lattice spacing is finite. It is 
certainly interesting to investigate if the Gromes relation is exactly satisfied 
in the continuum limit. Although we have investigated the gauge coupling at 
(3 = 6.3, we need further accuracy of the data at intermediate distances to 
achieve this. 

For V({r), if we only take into account the data for r > 0.25 fm at /3 = 6.0 
and, assuming that they are constant, we can fit them to a function 

K;fit = -^vl. (3.10) 

Due to the Gromes relation, we may identify this constant as a part of the 
string tension in V^{r). We then find a^^a^ = 0.0362(4) with Xmin/^df = 0.13, 
which is (77 ± 1) % of the string tension in VqIt). The corresponding fit curve 
is plotted in Fig. [HI 

While the Gromes relation is approximately satisfied, we find that the 
string tension a^i is not yet sufficient to reproduce the string tension crv2i- In 
other words there is still a missing amount of the string tension. We then notice 
that this must be supplied by V2(r). A fit of V2(r) to Eq. fl3.8p indeed leads 
to a^aa^ = 0.0070(7) and c = 0.288(7) with xLJ^di = 0.22. The fit curve is 
plotted in Fig. [91 Now the sum of a^i and crv2 reproduces (Tv2i- These findings 
suggest the existence of a long-ranged contribution in V2{r) whose magnitude 
is about one-fifth of a^i, which is (15 ± 2) % of the string tension in VqIt). 
We also attempted a fit with the expectation from perturbation theory, by 
simply fixing the string tension to be zero in the above fit. In this case the fit 
clearly fails, since Xmin/^df = 44. From the phenomenological point of view, 
one might prefer a simple parametrization like a = a^i and crv2 = |45j . 
but the results obtained here slightly differ from this expectation. We wish to 
point out, though, that V2{r) should further be investigated at distances larger 
than 0.7fm, in order to corroborate a non-vanishing value of V2(r) for r — > 00. 
We may note that Eq. (13.81) is not the only functional form for V-^ir). For 
instance, the function V2.f^^^{r) = c'/r^ with p = 1.51(4) and c'a^"^ = 0.205(9) 
also reproduces the data quite well, with Xmin/^df = 0.34. 

Next, we discuss the spin-spin potentials V^lr) and V4(r). We first examine 



° Here and in the following, we attach a subscript to a so as to distinguish the 
target function in the fit. 
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Vs^r) (see Fig. [T2]) if the ansatz motivated by one-gluon-exchange in Eq. (13.91) 
is appropriate. The fit to this function yields the coefficient c = 0.214(2) with 
Xmin/^df = 3.7. This value of Xmin/^df is relatively large and the result for c 
is 28 % smaller than the Coulombic coefficient in VqIt). A better fit can be 
achieved using an ansatz in which the power of 1/r is left as a free parameter, 
i.e. 



3r' 



(3.11) 
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This yields cV'^ = 0.171(10) and p = 2.80(6), where Xmin/^df = 0.79. The 
corresponding fit curve is plotted in Fig [121 The value of p is smaller than 3 
within 3 standard deviations. If one takes this result as face value, it indicates 
a deviation from the one-gluon-exchange potential. A deviation might actually 
be expected from the existence of the long-ranged contribution in V^'(r) and 
the relations in Table HI if we insert a function V{r) = —c/r + av2f into 
—V"{r) + V'{r)/r, we obtain 3c/r^ + (7^2/1^ at r 7^ 0. We have then examined 
if this function describes the data for given c and (Tv2, which are supplied 
from the V^'(r) fit. However, we have found that the resulting curve is not 
appropriate to describe the behavior of the data at all, since it lies above the 
data points at small distances. This tendency is practically due to the term 
1/r^, but this additional term 1/r also helps to lift the curve. It suggests that 
we need to add a negative contribution to such an ansatz. 

A possible candidate would then be a pseudo-scalar contribution, which 
is also closely related to the behavior of V^i^r) (see Fig. [T3]). In fact, if 
only the one-gluon-exchange interaction is considered in the vector kernel, 
2AV{r) = 2{V"{r) + 2V'{r)/r) leads to a 5-function as in Eq. ([33]), while 
if we insert the empirical behavior of V2{r), an additional term of 4(Tv2/t 
is generated for V4(r). Thus we expect a positive behavior at non-zero dis- 
tances. By contrast, the data is negative at small distances and almost zero 
for r > 0.2 fm. Let us now assume the presence of a pseudo-scalar interac- 
tion, P(r) = —g'e~^^^/r, where is the mass of the lightest pseudo-scalar 
particle, and g' is the corresponding effective coupling to quarks. This cer- 
tainly generates a negative contribution, AP(r) = —g'm?ge~^^^/r, to Vi{r). 
Note that the pseudo-scalar interaction P{r) is often used in the one-boson- 
exchange model for describing the nucleon-nucleon system, where pions play 
a relevant role [46]. In our simulation, however, since the effects of dynami- 
cal fermions are neglected due to our use of the quenched approximation, the 
lowest mass in the pseudo-scalar channel cannot be identified with the pion 
mass but rather with the lightest glueball mass. 

We have then performed a fit to 

VA-Ar) = -g'mf—-+A^, (3.12) 

where we have assumed mg = 2.47 GeV, which is taken from the recent lattice 
studies of the glueball masses [17], and treated g' and crv4 as free parameters. 
The resuh was g' = 0.292(12) and a^^a^ = 0.0015(3) with xLJ^di = 5.1, 
and the corresponding curve is put in Fig. [13] (if rrig is relaxed to be a free 
parameter, Xmin/^df is significantly reduced). We find that o"v4 (notice that 
this value is not zero) is not exactly a^2, but as the relation among interaction 
kernels is not exact but derived within the instantaneous approximation, such 
a deviation may occur, especially for the nonperturbative pieces. On the other 
hand, the value of g' is very close to c = 0.297(1) determined from the force 
(see Table ]3]). If we impose (Jv4 = 0"v2, we obtain here XminZ-^df ~ 100. However 
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Fig. 14. A possible structure of the spin-spin potential V3(r). The lattice data are 
the same as in Fig. [T21 The dotted (dashed) line corresponds to the vector (pseudo-s- 
calar) contribution, Eq. ()3.13p (Eq. ()3.14p ) and the sum of these two contributions 
is indicated by the solid line. 



. 00 




-0.25 



Fig. 15. r{2V3 + V4) - 6V^ at /? = 6.0. 



the apparent shape of the curve is not affected so much as is mostly governed 
by the glueball mass we set. 

Now we come back to V^lr) and examine if the sum of the (positive) vector 
and the (negative) pseudo-scalar contributions, with parameters estimated by 
the V4(r) fit, can describe the behavior of V3(r). Here, we consider the sum of 
two functions, 

Vr\r) = ^ + ^, (3.13) 
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Table 5 

Fit results of the spin-dependent potentials at (3 = 6.0 on the 20^40 lattice. (*) if 
we relax rUg to be a free parameter, Xmin/^^df is significantly reduced. 



Potential 


Fit ran^ 


;e (r/a) 


Fit function and parameters 


Xmin/^df 


Vlir) 






K = 






3 - 


- 7 


aa^ = 0.0362(4) 


0.13 








yfi',(r) = a + c/r2 






3 - 


- 6 


ao? = 0.0070(7), c = 0.288(7) 


0.22 








ykir) = d/r^ 






3 - 


- 6 


ddP-'^ = 0.205(9), p = 1.51(4) 


0.34 


Vi{r)-Vl{r) 






Vi,{r)=cj + c/r^ 






3 - 


- 7 


oo? = 0.0426(9), c = 0.293(9) 


0.26 








yfit(r) = 3c/r3 






3 - 


- 7 


c = 0.214(2) 


3.7 








^fit(r) = 3c7rP 






3 - 


- 7 


c'aP-^ = 0.171(10), p = 2.80(6) 


0.79 








Vfi_t{r) = -g'm?ge~'^a'^ /r + Aa/r 






2 - 


- 7 


g' = 0.292(12), aa^ = 0.0015(3), 
rUga = 1.16 (fixed) 


5.1* 



V<«V) = P"ir) - ^ = -/(4 + ^ + '^)e-' . (3.14) 

where c is taken from V(^{r). It is interesting to see Fig. [T3] that the curve 
^3*'^''('") + '^3^^\^)^ which is plotted with the solid line, can go through the 
data at P = 6.0. 

In most of previous works, it was concluded that V^lr) is consistent with 
a ^-function, which may only be true at very small distances. However, as 
we demonstrated, the behavior of Vs{r) and V4(r) at distances r > 0.2 fm 
can be consistently explained by assuming the existence of the pseudo-scalar 
contribution as well as the vector contribution. Note furthermore that the 
combination of the potentials r(2V3 + V4) — 6V2 should be zero at non-zero 
distances, if the interaction kernel contains only the pure vector component 
without a linear term and no pseudo-scalar contribution (see Table Hj) [23] . 
However, as shown in Fig. [15], we find that this combination is non-vanishing 
within our accuracy, so that some of these assumptions are probably not ap- 
plicable. Of course, our discussion on the pseudo-scalar contribution is as yet 
speculation, which needs to be checked in future works. A possible way of 
doing this is to investigate V4(r) in the presence of dynamical quarks (pions) 
in full QCD simulations, and to examine whether one can indeed observe a 
behavior like cx —e~"^'^^/r for sufficiently small quark masses. Some of previ- 
ous works in Refs. [271128] have been carried out in full QCD, but the data 
quality is not sufficient to draw any conclusion. In any case, we expect to raise 
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further discussions on the structure of the spin-spin potentials. 

To close the discussion on the functional form, we note that the Gromes 
inequalities V^^r) > V^^r) and 2V3(r) + V4(r) > are certainly satisfied. For 
instance, the latter inequality is immediately checked through 2Vs{r) + V4{r) = 
6 /q°° dt{{g^Bx(0, 0)Brc{f, t))), which is positive at all available r as can be seen 
from Tables IA.3tlA.5l in the Appendix. We summarize all fit results of the 
functional form in Table [51 

4 Summary 

We have investigated the spin- dependent corrections to the static po- 
tential at 0(l/m^) in SU(3) lattice gauge theory. These corrections, usually 
called the spin-dependent potentials, are represented as the integral of the 
field strength correlators on the quark- ant iquark source with respect to the 
relative temporal distance between two field strength operators. We have used 
the Polyakov loop correlation function as the quark-antiquark source, and by 
employing the multi-level algorithm, we have obtained remarkably clean data 
for the expectation values of the field strength correlators and, in turn, for the 
spin-dependent potentials up to intermediate distances of around r ~ 0.6 fm. 
The spectral representation of the field strength correlator in a finite peri- 
odic volume has been exploited in order to extract the potential with less 
systematic error. 

The observation we have made for the spin-dependent potentials in 
Eq. fll.ip is as follows. The spin-orbit potential V({r) is clearly long-ranged, is 
negative at all distances and constant at r > 0.25 fm. The other spin-orbit po- 
tential V2 (r) is positive at all distances and shows a behavior decreasing with r. 
However, it has a finite tail up to intermediate distances, which cannot be ex- 
plained at least by the one-gluon-exchange interaction. The Gromes relation 
V^'(r) = V2(r) — V({r) is satisfied within (8 ± 1) % accuracy at intermediate 
distances in the present simulation. Within this relation, the constant value 
in V({r) reproduces (77 ± 1) % of the string tension in lo'(r) and (15 ± 2) % 
of the string tension are found to be supphed by ^'(r). The spin-spin (ten- 
sor) potential V^i^r) is positive at all distances and is decreasing as a function 
of r. The behavior is slightly more moderate than the expectation of the one- 
gluon-exchange picture oc 1/r'^. The other spin-spin potential V^lr) exhibits 
a negative short-ranged behavior. This short-ranged behavior, as well as the 
behavior of V^^r), could be explained if the exchange of the pseudo-scalar 
glueball is assumed in addition to the one-gluon-exchange type interaction. 

In this paper we have not carried out a detailed comparison of the lattice 
result of the spin-dependent potentials with perturbation theory, e.g. along 
the lines of Necco and Sommer for the static potential [37^8] . In this sense, 
although we have observed a certain deviation from the expectation of leading 
order perturbation theory at intermediate distances, it is not yet clear that 
from which distance a perturbative description becomes inadequate. Clearly, 
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it requires further systematic studies, where the renormahzation of the field 
strength operator and also the matching coefficients are worth to be reconsid- 
ered. However, we expect that the numerical procedures we have demonstrated 
in this paper is quite useful for such a work. Then, it is interesting to use the 
result as inputs of phenomenological models P9|50|51] or to compare with the 
various QCD vacuum models [52|53II54] . It may be interesting to note that 
the existence of a long-ranged contribution in V2{r) is suggested in Ref. |54j . 
independently of the present work. 

Finally we note that our numerical procedures are also applicable to the 
evaluation of other relativistic corrections like the velocity-dependent poten- 
tials [55|56llllj and the potential at 0(l/m) [14fl5llllll57f58ll59ll60f61j . which 
are also represented as the field strength correlators on the quark- ant iquark 
source with different combination of the field strength operators from the spin- 
dependent potentials. The first lattice result on the potential at 0(l/m) was 
published in Ref. [M] . 
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A Collection of numerical values 

A.l Tree-level improvement of the quark- antiquark distances 
Table A.l 

The quark- antiquark distances for the static potential Vo{rj), the force lo'(f) (or 
VQ'(fc)) and the second derivative V^"(r) with the tree- level improvement |38|37|32] . 



r/a 


ri/a 


r/a 


rc/a 


r/a 


1 


0.925 








2 


1.855 


1.358 


1.649 


1.788 


3 


2.889 


2.277 


2.654 


2.700 


4 


3.922 


3.312 


3.729 


3.729 


5 


4.942 


4.359 


4.794 


4.786 


6 


5.954 


5.393 


5.837 


5.833 


7 


6.962 


6.414 


6.865 


6.864 


8 


7.967 


7.428 


7.885 


7.886 


9 


8.971 


8.438 


8.899 


8.901 



A. 2 HM factors 
Table A.2 

The HM renormalization factors at /? = 6.0 on the 20'^ lattice (upper) and {3 = 6.3 
on the 24'^ lattice (lower), where the quark-antiquark system is set along the x axis. 
Thus, one should observe = Ze^-, = Zb^- 



r/a 




ZEy 






ZBy 




2 


1.59446(4) 


1.63031(6) 


1.63038(5) 


1.67833(16) 


1.67614(12) 


1.67600(15) 


3 


1.61170(4) 


1.62498(6) 


1.62503(5) 


1.67764(16) 


1.67661(12) 


1.67651(15) 


4 


1.61620(4) 


1.62338(6) 


1.62339(6) 


1.67735(16) 


1.67676 (12) 


1.67669(14) 


5 


1.61777(6) 


1.62282(6) 


1.62282(7) 


1.67726(16) 


1.67687(13) 


1.67678(15) 


6 


1.61846(6) 


1.62250(7) 


1.62262(6) 


1.67721(16) 


1.67695(13) 


1.67683(16) 


7 


1.61877(10) 


1.62246(8) 


1.62233(8) 


1.67726(16) 


1.67684(13) 


1.67680(15) 


8 


1.61879(15) 


1.62225(17) 


1.62232(16) 


1.67708(18) 


1.67682(19) 


1.67674(19) 


2 


1.54232(3) 


1.56529(3) 


1.56526(3) 


1.60307(7) 


1.60179(9) 


1.60185(7) 


3 


1.55417(5) 


1.56151(4) 


1.56154(4) 


1.60271(7) 


1.60217(10) 


1.60226(7) 


4 


1.55717(6) 


1.56048(4) 


1.56049(6) 


1.60266(7) 


1.60225(10) 


1.60235(7) 


5 


1.55810(9) 


1.56013(6) 


1.56004(7) 


1.60260(9) 


1.60231(11) 


1.60231(10) 


6 


1.55853(10) 


1.55974(11) 


1.55986(13) 


1.60247(11) 


1.60229(12) 


1.60253(11) 


7 


1.55829(18) 


1.55921(23) 


1.55980(14) 


1.60238(18) 


1.60216(17) 


1.60223(16) 


8 


1.55855(40) 


1.55974(37) 


1.55975(32) 


1.60260(29) 


1.60249(28) 


1.60217(23) 


9 


1.55961(43) 


1.55954(64) 


1.56021(48) 


1.60324(36) 


1.60324(39) 


1.60331(44) 
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A. 3 Fit results of the field strength correlators 



Table A.3 

Fit results of the field strength correlators with Eqs. (j2.8p - ()2.10p at /? = 6.0 on the 
20"^ lattice, n-max is the maximum truncation level of the spectral representation. 



C{t) 


r/a 


Fit rang 


;e (t/a) 




a2 /~ dt tC{t) 


Xmin/^df 




2 


2 - 


10 


2 


-0.01801(4) 


0.43 




3 


2 - 


10 


2 


-0.01808(6) 


0.99 




4 


2 - 


10 


2 


-0.01805(10) 


0.30 




5 


2 - 


10 


2 


-0.01800(19) 


0.45 




6 


2 - 


10 


2 


-0.01795(32) 


1.9 




2 


1 - 


10 


3 


0.03618(3) 


2.0 




3 


1 - 


10 


3 


0.01948(4) 


2.7 




4 


1 - 


10 


3 


0.01265(10) 


0.81 




5 


1 - 


10 


2 


0.00917(14) 


2.0 




6 


1 - 


10 


2 


0.00795(42) 


0.88 


C{t) 


r/a 


Fit rang 


;e (t/a) 




a3 /~ dt C{t) 




{{g^B^{0,0)B^{r,t))) 


2 


1 - 


10 


3 


0.01648(3) 


2.2 




3 


1 - 


10 


3 


0.00743(3) 


0.97 




4 


1 - 


10 


3 


0.00381(3) 


0.85 




5 


1 - 


10 


2 


0.00217(3) 


2.0 




6 


1 - 


10 


2 


0.00128(5) 


1.8 


{{g^By{0,0)By{r,t))) 


2 


1 - 


10 


3 


-0.01226(3) 


0.77 




3 


1 - 


10 


3 


-0.00443(3) 


0.40 




4 


1 - 


10 


3 


-0.00155(3) 


0.37 




5 


1 - 


10 


3 


-0.00067(3) 


0.46 




6 


1 - 


10 


3 


-0.00033(4) 


0.48 
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Tabic A.4 

Fit results of the fiekl strength eorrelators at 6 — 6.0 on the 20'^ 40 lattiee. 



C(t) 


via 


Fit range it la) 




f~ tCit) 


/vmin/ ui 


\\y yv^^ ^ j-^zy^") ^) II 


2 


2-20 


2 


-0.01798(13) 


2.7 




3 


2-20 


2 


-0.01809(20) 


2.5 




4 


2-20 


2 


-0.01800(27) 


2.9 




5 


2-20 


2 


-0.01809(25) 


1.4 




6 


2-20 


2 


-0.01831(40) 


1.8 




7 


2-20 


2 


-0.01853(137) 


2.6 


((q'^BJO, 0)EJr, t))) 
\\u yy^^ ^}^z\' 5 ^) II 


2 


1-20 


3 


0.03618(5) 


1.4 




3 


1-20 


3 


0.01950(12) 


5.0 




4 


1-20 


2 


0.01253(18) 


1.9 




5 


1-20 


2 


0.00920(31) 


2.4 




6 


1-20 


2 


0.00703(80) 


5.7 




7 


1-20 


2 


0.00457(68) 


2.6 


C{t) 


r/a 


Fit range {tja) 






Xmin/^df 


{{g^B,{0,0)B,ir,t))) 


2 


1-20 


3 


0.01642(7) 


1.5 




3 


1-20 


3 


0.00742(5) 


2.2 




4 


1-20 


3 


0.00372(7) 


1.5 




5 


1-20 


3 


0.00206(6) 


1.5 




6 


1-20 


2 


0.00133(16) 


5.5 




7 


1-20 


2 


0.00071(18) 


4.5 


{{9^By{0,0)By{r,t))) 


2 


1-20 


3 


-0.01226(6) 


3.1 




3 


1-20 


3 


-0.00442(7) 


2.3 




4 


1-20 


3 


-0.00166(7) 


2.8 




5 


1-20 


3 


-0.00069(6) 


2.6 




6 


1-20 


3 


-0.00029(9) 


1.8 




7 


1-20 


3 


-0.00030(39) 


1.9 
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Table A.5 

Fit results of the field strength eorrelators at /3 = 6.3 on the 24^ lattice. (*)one of 
the excitation energies in the expansion was fixed so as to make the fit stable. 



C{t) 


r/a 


Fit ranj 


;e {t/a) 




tC(i) 




{{g^By{0,0)E,{0,t))) 


2 


2 - 


12 


2 


-0.00980(10) 


1.1 




3 


2 - 


12 


2 


-0.00872(29) 


1.2 




4 


2 - 


12 


2 


-0.00768(30) 


1.4 




5 


2 - 


12 


2 


-0.00727(28) 


4.4* 




6 


2 - 


12 


2 


-0.00664(110) 


0.94 


{{g^ByiO,0)E,{r,t))) 


2 


1 - 


12 


3 


0.03145(7) 


0.28 




3 


1 - 


12 


2 


0.01620(16) 


3.0 




4 


1 - 


12 


2 


0.01008(37) 


2.0 




5 


1 - 


12 


2 


0.00632(43) 


0.85 




6 


1 - 


12 


2 


0.00446(109) 


2.3 


C{t) 


r/a 


Fit ranj 


?e {t/a) 






Xmin/-^df 




2 


1 - 


12 


3 


0.01452(3) 


0.49 




3 


f - 


12 


3 


0.00653(3) 


4.4 




4 


1 - 


12 


2 


0.00332(6) 


1.2 




5 


1 - 


12 


2 


0.00188(8) 


0.87 




6 


1 - 


12 


1 


0.00123(12) 


0.49 


{{g'By{0,0)By{r,t))) 


2 


1 - 


12 


3 


-0.01203(3) 


0.95 




3 


1 - 


12 


3 


-0.00436(3) 


1.4 




4 


1 - 


12 


3 


-0.00168(7) 


1.5 




5 


1 - 


12 


2 


-0.00084(18) 


4.8 




6 


1 - 


12 


1 


-0.00046(12) 


1.8 
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A. 4 Spin- dependent potentials 



Table A.6 

The spin-dependent potentials at /? = 6.0 on the 20"^ lattice. 



r/a 








a^V4 


2 


-0.03603(8) 


0.07235(6) 


0.05749(8) 


-0.01607(13) 


3 


-0.03616(12) 


0.03896(9) 


0.02373(8) 


-0.00286(14) 


4 


-0.03609(20) 


0.02531(19) 


0.01072(8) 


0.00143(17) 


5 


-0.03600(39) 


0.01836(28) 


0.00568(9) 


0.00165(14) 


6 


-0.03591(65) 


0.01590(84) 


0.00322(13) 


0.00123(20) 


Table A.7 










The spin-dependent potentials at /3 = 6.0 on the 20^40 lattice. 


r/a 










2 


-0.03596(25) 


0.07235(11) 


0.05736(19) 


-0.01620(27) 


3 


-0.03618(41) 


0.03900(25) 


0.02368(20) 


-0.00284(29) 


4 


-0.03600(54) 


0.02507(36) 


0.01075(21) 


0.00081(30) 


5 


-0.03619(50) 


0.01840(61) 


0.00549(18) 


0.00138(27) 


6 


-0.03662(81) 


0.01405(160) 


0.00323(37) 


0.00152(47) 


7 


-0.03706(274) 


0.00914(136) 


0.00203(103) 


0.00021(137) 


Table A.8 










The spin-dependent potentials at /? = 6.3 on the 24"^ lattice. 




r/a 










2 


-0.01960(19) 


0.06290(15) 


0.05310(9) 


-0.01910(15) 


3 


-0.01744(58) 


0.03240(31) 


0.02178(9) 


-0.00439(14) 


4 


-0.01539(61) 


0.02017(74) 


0.01000(16) 


-0.00009(31) 


5 


-0.01443(79) 


0.01265(87) 


0.00546(41) 


0.00038(70) 


6 


-0.01329(221) 


0.00893(217) 


0.00339(37) 


0.00057(52) 



33 



References 

[1] W. Lucha, F.F. Schoberl and D. Gromes, Phys. Kept. 200 (1991) 127. 

[2] W. Buchmiiller (Ed.), Quarkonia, Current physics-sources and comments, 
vol. 9, (North-Holland, 1992). 

[3] G.S. Bah, Phys. Kept. 343 (2001) 1, hep-ph/0001312j 

[4] N. BrambUla et al., Heavy quarkonium physics, CERN YeUow Report (2005), 



hep-ph/0412158 



[5] E. Eichten and F. Feinberg, Phys. Rev. Lett. 43 (1979) 1205. 

[6] E. Eichten and F. Feinberg, Phys. Rev. D23 (1981) 2724. 

[7] M.E. Peskin, Aspects of the dynamics of heavy quark systems, 11th Int. SLAG 
Summer Inst, on Particle Physics, Stanford, CA, Jul 18-26, 1983, SLAC-PUB- 
3273. 

[8] D. Gromes, Z. Phys. G22 (1984) 265. 

[9] Y.J. Ng, J.T. Pantaleone and S.H.H. Tye, Phys. Rev. Lett. 55 (1985) 916. 



[lo; 

[11 
[12 

[13: 

[14 
[15 

[16; 
[ir 

[18 
[19 

[2o; 

[21 
[22; 



J.T. Pantaleone, S.H.H. Tye and Y.J. Ng, Phys. Rev. D33 (1986) 777. 



A. Pineda and A. Vairo, Phys. Rev. D63 (2001) 054007, |hep-ph70009145 
Eiiatum-ibid D64 (2001) 039902. 

W.E. Gasweh and G.P. Lepage, Phys. Lett. B167 (1986) 437. 

G.T. Bodwin, E. Braaten and G.P. Lepage, Phys. Rev. D51 (1995) 1125, 



|hep-ph/ 9407339 1 Eiiatum-ibid D55 (1997) 5853. 

N. Brambilla, A. Pineda, J. Soto and A. Vairo, Nucl. Phys. B566 (2000) 275, 
|hep-ph/9907240j 



N. Brambilla, A. Pineda, J. Soto and A. Vairo, Phys. Rev. D63 (2001) 014023, 
[hep^ph7 0002250i 

N. Brambilla, A. Pineda, J. Soto and A. Vairo, Rev. Mod. Phys. 77 (2005) 



1423, |h ep-ph/0410047 



A.V. Manohar, Phys. Rev. D56 (1997) 230, |hep-ph/9701294 
D. Gromes, Z. Phys. G26 (1984) 401. 

N. Brambiha, D. Gromes and A. Vairo, Phys. Rev. D64 (2001) 076010, 



hep-ph/0104068, 



Ph. de Forcrand and J.D. Stack, Phys. Rev. Lett. 55 (1985) 1254. 
G. Michael and P.E.L. Rakow, Nucl. Phys. B256 (1985) 640. 
C. Michael, Phys. Rev. Lett. 56 (1986) 1219. 



34 



[23; 

[24 

[25 

[26; 
[27; 
[28; 

[29 

[3o; 

[31 
[32 
[33; 
[34; 

[35; 

[36; 

[37; 
[38; 
[39; 
[4o; 

[41 
[42 
[43; 
[44 
[45 
[46; 
[47 
[48; 



A. Huntley and C. Michael, Nucl. Phys. B286 (1987) 211. 

M. Campostrini, K. Moriarty and C. Rebbi, Phys. Rev. Lett. 57 (1986) 44. 

M. Campostrini, K. Moriarty and C. Rebbi, Phys. Rev. D36 (1987) 3450. 

I.J. Ford, J. Phys. G15 (1989) 1571. 

Y. Koike, Phys. Lett. B216 (1989) 184. 

K.D. Born, E. Laermann, T.F. Walsh and P.M. Zerwas, Phys. Lett. B329 (1994) 
332. 

G.S. Bah, K. Schilling and A. Wachter, Phys. Rev. D55 (1997) 5309, 
hep-lat/9611 025| 



G.S. Bah, K. Schilling and A. Wachter, Phys. Rev. D56 (1997) 2566, 
[hep4a:t79703019| 



M. Liischer and P. Weisz, JHEP 09 (2001) 010, hep-lat/0108014 



M. Liischer and P. Weisz, JHEP 07 (2002) 049, hep-lat/0207003 



M. Koma, Y. Koma and H. Wittig, PoS LAT2005 (2005) 216, |hep-lat /0510059[ 



Y. Koma, M. Koma and H. Wittig, Phys. Rev. Lett. 97 (2006) 122003, 



hep-lat/0607009 



G.P. Lepage and P.B. Mackenzie, Phys. Rev. D48 (1993) 2250,|hep-lat/9209022| 

Y. Koma, M. Koma and P. Majumdar, Nucl. Phys. B692 (2004) 209, 
|hep-lat/03110T6| 



S. Necco and R. Sommer, Nucl. Phys. B622 (2002) 328, |hep-lat70108008f 



R. Sommer, Nucl. Phys. B411 (1994) 839, |hep-lat/9310022 

M. Liischer, K. Symanzik and P. Weisz, Nucl. Phys. B173 (1980) 365. 

M. Liischer, Nucl. Phys. B180 (1981) 317. 

N.D. Hari Dass and P. Majumdar, JHEP 10 (2006) 020, hep-lat/0608024] 
D. Gromes, Phys. Lett. B202 (1988) 262. 
D. Gromes, Nucl. Phys. B131 (1977) 80. 

F. Gesztesy, H. Grosse and B. Thaller, Phys. Rev. D30 (1984) 2189. 
W. Buchmiiller, Phys. Lett. B112 (1982) 479. 

F. Gross, J.W. Van Orden and K. Holinde, Phys. Rev. C45 (1992) 2094. 
Y. Chen et al., Phys. Rev. D73 (2006) 014516, |hep4at7o'5l0074t 



S. Necco and R. Sommer, Phys. Lett. B523 (2001) 135, hep-ph/0109093 



35 



[49; 
[5o; 

[51 

[52 
[53; 

[54; 
[55 

[56; 
[57; 

[58; 
[59; 

[60; 

[61 



D. Ebert, V.O. Galkin and R.N. Faustov, Phys. Rev. D57 (1998) 5663, 
|hep-ph/9712318[ 

D. Ebert, R.N. Faustov and V.O. Galkin, Phys. Rev. D62 (2000) 034014, 
|hep-ph/9911283| 

D. Ebert, R.N. Faustov and V.O. Galkin, Phys. Rev. D67 (2003) 014027, 



hep-ph/0210381[ 



Y.A. Simonov, Nucl. Phys. B324 (1989) 67. 



N. Brambilla and A. Vairo, Phys. Rev. D55 (1997) 3974, hep-ph/9606344 



F. Jugeau and H. Sazdjian, Nucl. Phys. B670 (2003) 221, hep-ph/0305021 



A. Barchielh, E. Montaldi and G.M. Prosper!, Nucl. Phys. B296 (1988) 625, 
Erratum-i6id B303 (1988) 752. 

A. Barchielh, N. Brambilla and G.M. Prosperi, Nuovo Cim. A103 (1990) 59. 

K. Melnikov and A. Yelkhovsky, Nucl. Phys. B528 (1998) 59, |hep-ph/ 9802379| 

A. H. Hoang, Phys. Rev. D59 (1999) 014039, |hep-ph/9803454| 

N. Brambilla, A. Pineda, J. Soto, and A. Vairo, Phys. Lett. B470 (1999) 215, 
|hep-ph/9910238[ 



B.A. Kniehl, A. A. Penin, M. Steinhauser and V.A. Smirnov, Phys. Rev. D65 
(2002) 091503, |hep-ph/0106135[ 

B.A. Kniehl, A. A. Penin, V.A. Smirnov and M. Steinhauser, Nucl. Phys. B635 



(2002) 357, hep-ph/0203166 



36 



